Matrix product states for critical spin chains: 
finite size scaling versus finite entanglement scaling 
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We investigate the use of matrix product states (MPS) to approximate ground states of critical 
quantum spin chains with periodic boundary conditions (PBC). We identify two regimes in the 
(N, D) parameter plane, where N is the size of the spin chain and D is the dimension of the 
MPS matrices. In the first regime MPS can be used to perform finite size scaling (FSS). In the 
complementary regime the MPS simulations show instead the clear signature of finite entanglement 
scaling (FES). In the thermodynamic limit (or large TV limit), only MPS in the FSS regime maintain 
a finite overlap with the exact ground state. This observation has implications on how to correctly 
perform FSS with MPS, as well as on the performance of recent MPS algorithms for systems with 
PBC. It also gives clear evidence that critical models can actually be simulated very well with MPS 
by using the right scaling relations; in the appendix, we give an alternative derivation of the result 
of Pollmann et al. [Phys. Rev. Lett. 102, 255701 (2009)] relating the bond dimension of the MPS 
to an effective correlation length. 

PACS numbers: 02.70.-c, 03.67.-a, 05.10.Cc, 75.10.Pq 



I. INTRODUCTION 

Quantum many body systems are very hard to study 
due to the exponential growth of their Hilbert space with 
the number of constituents. One possible cure to this 
issue for one dimensional systems is to describe their 
ground states as matrix product states (MPS) HH2|. This 
family of states is known to be well suited to study gaped 
ID phases [4j where for generic systems almost exact re- 
sults can be obtained with matrices whose size does not 
depend on the size of the system. Even more, for sev- 
eral gapped ID systems the exact ground state can be 
expressed in terms of translationally invariant MPS with 
very small bond dimension Gapless ID phases are 

harder to simulate with MPS since the size of the ma- 
trices necessary to obtain good approximations of their 
ground states increases polynomially with the size of the 
system. This is particularly unfortunate since the univer- 
sal low energy information encoded in the gapless phase 
becomes apparent only for large systems. 

Luckily such universal information is also encoded in 
the way a state approaches the thermodynamic limit and 
one can extract it by using the celebrated finite size scal- 
ing (FSS) technique i6i]. This technique amounts to study 
larger and larger systems in a gapless phase (that due to 
the finite size of the system becomes gapped) and ex- 
tract universal properties through the dependence of the 
observables on the system size. 

In the context of MPS, one can use an alternative ap- 
proach to study gapless phases. It is called finite entan- 
glement scaling (FES) [7] and amounts to study the scal- 
ing of the expectation value of observables in the ground 
state of infinite chains described by MPS with fixed bond 
dimension and thus finite entanglement |S]- 

Both the existence of FSS and FES close to a con- 



formal fixed point are a direct consequence of conformal 
invariance [9l [10] . If IV is the chain length and D the 
MPS bond dimension, then FSS corresponds to taking 
D — » oo first and then taking N — > oo, whereas FES 
consists in taking N — > oo first and then D — > oo. 

An important question to ask is whether FSS and FES 
provide the same universal information. Since the pro- 
posal of FES for simulations with MPS [7] it has been 
shown that indeed quantities such as critical exponents 
related to local observables or the central charge of the 
model can be extracted with the help of this technique 
|llH17j . in a similar way as it is normally done with FSS 
techniques. Here we will show, however, that some care 
is required in order to differentiate between the effects of 
FES and those of FSS. 

Specifically, we consider critical systems with periodic 
boundary conditions (PBC), and describe their transla- 
tionally invariant ground states using translationally in- 
variant MPS. In order to perform FSS one should obtain 
for each system size iV a sequence of increasingly ac- 
curate MPS approximations with growing bond dimen- 
sion D, which for large enough D converges to the exact 
ground state. Importantly, we find that for an intermedi- 
ate range of values of D, for which local observables are 
already reproduced with high accuracy and show scal- 
ing behavior, the MPS approximation is almost com- 
pletely orthogonal to the exact ground state (resulting 
e.g. in failure to reproduce correlation functions at dis- 
tance N/2, as previously illustrated in the inset of Fig. 
11 and 12 of Ref. [E]). In other words, reasonably 
converged values of (and/or scaling behavior for) local 
observables including the ground state energy, are not 
sufficient criteria to establish that some MPS is a good 
approximation to the ground state of a critical PBC sys- 
tem. Instead, in order to properly apply FSS, for each 
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system size N one should consider MPS with a bond di- 
mension D larger than some threshold value -Do, where 
D depends both on N and on the spin model. 

Our results have important consequences for the de- 
sign of algorithms that simulate PBC chains with MPS. 
Simulating PBC systems with MPS is computationally 
much more expensive |19j than simulating the same sys- 
tem with open boundary conditions (OBC). Nonetheless, 
when studying critical ground states, systems with PBC 
are known to approach the thermodynamic limit much 
faster than systems with OBC, and therefore they offer 
a much better framework for FSS. For this reason, sub- 
stantial effort l2"0Tl2"2"] has been made to try to lower 
the computational cost of MPS simulations with PBC. 
Two types of approaches have been pursued. One con- 
sists in building a MPS for a finite system with PBC by 
using the translationally invariant MPS tensor that has 
been optimized in an infinite chain with OBC |21j . This 
approach is equivalent to a crude approximation of the 
MPS transfer matrix: the D 2 x D 2 matrix is approxi- 
mated only by its dominant eigenvector |23j . The second 
approach 18 , 20, 22 accounts for PBC by retaining more 
than one eigenvector in the approximation of the trans- 
fer matrix. We show in this work that the first approach 
fails to provide an accurate ground state approximation 
for critical PBC systems. A detailed comparison of these 
algorithms can be found in Appendix [C] 

We will build our arguments by studying two paradig- 
matic critical spin chains: the quantum Ising model (IS) 
and the quantum Heisenberg (HB) model, for chains with 
PBC and linear size N. The ground states are encoded 
in MPS of a given bond dimension D. Even if the Hamil- 
tonian is critical, both the finite size of the chain and the 
finite bond dimension of the MPS induce a gap 

= & = ^ (1) 



AE D = & cx D- K (2) 

where x\ is the smallest critical exponent of the theory 
[Mj and k is the exponent for the scaling of the effective 
correlation length of MPS simulations with finite bond 
dimension [Jj [10] . Depending on which of the two gaps 
dominates, the system is in one of the two regimes 

£,d > &v ■ FSS regime (3) 



6v > : FES regime (4) 

The presence of two regimes in the PBC chain can be 
intuitively understood in the following way: in the FES 
regime defined by equation Q the small dimension of the 
MPS matrices implies that the system is not aware of its 
geometry. Thus the boundaries do not play any role. 



In the FSS regime, defined by equation ^ on the other 
hand, the size of the matrices is big enough to notice the 
presence of the boundaries and thus different choices of 
boundary conditions lead to different MPS. 
For simulations where 

6v =s Zd (5) 

we find for all values of N and D that are accessible nu- 
merically the presence of an abrupt transition between 
the FSS and FES regimes (for a related work see also 
Ref. [25]). One way to observe this transition is by 
looking at the difference between the exact ground state 
energy in the thermodynamic limit and the energy of 
MPS approximations with different N and D. For fixed 
D these plots show a steep transition between the FSS 
regime where the difference scales like cx N~ 2 to the 
FES regime where the difference does not depend on 
N. Another way is to look at the overlap between MPS 
with different D for fixed chain length TV: starting off 
with a MPS with some big D max , we look at its overlap 
with MPS with decreasing D. We then observe how the 
initially smoothly decreasing overlap abruptly drops to- 
wards lower values at some D r , unambiguously showing 
the transition to the FES regime. Now the overlap is a 
global variable and as such indeed aware of the boundary 
conditions. The main finding of this paper is that states 
in the FES regime, while possessing the same local uni- 
versal properties [Jj like those in the FSS regime, turn 
out to have vanishing overlap with them. 

We also present a possible technique to determine if 
a given bond dimension is sufficient to enter the FSS 
regime, so that we can give the computationally most fa- 
vorable recipe to access global universal properties that 
depend on the boundaries (for a discussion of these prop- 
erties see e.g. [26]). 

The paper is organized as follows. We start by intro- 
ducing the IS and HB models as well as the technique 
used to simulate them in section HU In section III Al we 
present numerical evidence for the presence of the FSS 
and FES regimes in MPS simulations of PBC chains by 
looking at the ground state energy. In section [TlB| we dis- 
cuss how to identify the sharp transition between the two 
regimes by looking at the overlap. Then, in section al B 1[ 
we give a more detailed view on the transition between 
the two scaling regimes. Section |II C| gives a recipe for 
obtaining the minimal bond dimension needed to observe 
global universal properties of critical systems. In section 
|IID| we perform a numerical study of the transition and 
for the IS model we can provide evidence for its persis- 
tence in the thermodynamic limit. For the HB model 
we are not able to do the same due to the coarser preci- 
sion of our simulations for this model. In section Hi El we 
provide a numerical analysis of the scaling function for 
the energy difference that reveals a two-parameter scal- 
ing analogous to the one found in the context of critical 
2D classical systems by Nishino in Ref. [27] . We conclude 
with a discussion of the implications of our results and 
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with a brief outline of future developments. 

All technical details are contained in the appendices. 
There we first provide an alternative way to derive the 
analytical result for the scaling exponent k in £d oc D k , 
which we find more intuitive than the one given in [TO] . 
Then we show how our algorithm can be used in order to 
extract k from the numerical results for the ground state 
energy. For the IS model, we are able to provide a nu- 
merical confirmation for the persistence of the transition 
between the FSS and the FES regime in the thermody- 
namical limit. 



II. NUMERICAL RESULTS 

We will use throughout this work the algorithm pre- 
sented in [18]. That algorithm exploits the translational 
invariance of the models we study by using an ansatz 
based on translationally invariant MPS. This means that 
the MPS tensors at each site of the chain are identical 
thus reducing the cost of the simulation by a factor N. 
The energy is minimized by means of a conjugate gra- 
dient method in the subspace spanned by real and sym- 
metric MPS with bond dimension D. The computational 
cost scales like 0(mnD 3 ) + 0(n 2 D 3 ). m and n are pa- 
rameters whose magnitude depends on the entanglement 
of the ground state of the model under consideration. For 
more details on the method we refer the reader to that 
work. The two paradigmatic models we have considered 
are the critical quantum Ising model described by the 
Hamiltonian 



N N 
i=l i=l 

and the Heisenberg model described by the Hamiltonian 

N N 
i=l i=l 

(7) 

where the ID lattice is considered to be periodic. Both 
Hamiltonians are critical which means that their gap be- 
tween the ground state and the first excited state closes 
as an inverse power of N as described in Eq. [I] 

As a matter of fact the analysis in this work was trig- 
gered by a comprehensive study of the precision of the 
algorithm presented in [TH]. We originally wanted to as- 
sess the usability of that method and to this end we simu- 
lated a plethora of different configurations {N, D} for the 
IS and the HB models. A selection of these simulation 
results is shown in Figure [T] where we plot the relative en- 
ergy precision of simulations with different D for many 
different chain lengths N. 

The shape of curves with constant D is very surpris- 
ing since it shows a fundamental deviation from what we 
would have expected. Our expectation was that for short 
chains the precision will be generally better than for long 
chains and that as N gets bigger and bigger, the precision 
will eventually saturate from below to the value obtained 
with the corresponding D when simulating the chain in 
the thermodynamic limit. Obviously the small N and 
the big N regimes are in accordance with our expecta- 
tion. However at some point between these limits we see 
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the emergence of a hump which indicates that something 
interesting is happening in that region. As a matter of 
fact we can show that if we interpret the position of the 
hump as an indicator for the effective correlation length 
£d of MPS with finite D, we can reproduce the theoreti- 
cally predicted result for the scaling of £d EE] with very 
good accuracy (see Appendix [A| . The results presented 
in this work provide the general framework to understand 
the emergence of the hump and to explain what is hap- 
pening when moving from the left to the right side of Fig. 

m 



A. Two different regimes for MPS simulations 

As already mentioned in the introduction a MPS sim- 
ulation close to the critical point is an example of a 
two scale problem. This is not something unexpected 
as it has been already pointed out in the context of 2D 
classical systems studied with the corner transfer ma- 
trix by Nishino and coworkers [57] and in the context of 
quantum phase transitions in ID quantum chains with 
open boundaries by one of the authors (section IIIG of 
Ref. [7J). In the scenario we are considering, the two 
scales appearing are i) the correlation length induced by 
the finite size of the system N of Eq. [I] and ii) the cor- 
relation length induced by the size of the matrices D of 
Eq. [2] Depending on the relation among the two stated 
in Eqs. [3] and [4] the system will be in one of the two dif- 
ferent regimes, respectively the FSS regime or the FES 
regime. 

The approach to the thermodynamic limit of the 
ground state energy, as function of the relevant parame- 
ter N or D, is very different in the two regimes so that 
we can use it as a footprint for them. In the FSS regime 
indeed it obeys the celebrated result by Cardy and Af- 
fleck [SHI from conformal field theory (CFT), 
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6N 2 



(8) 



where Eo(oo) is the thermodynamic limit, Vf is the Fermi 
velocity and c is the central charge of the considered 
model. In the thermodynamic limit, several authors have 
reported that |?1ITCHI2"5] 




FIG. 2. (Color online), a) Quantum Ising model, two 
regimes for simulations of PBC chains with lengths in the 
range 20 < N < 4 ■ 10 4 , and D in the range 4 < D < 64. 
Each line represents simulations performed at fixed D and 
different N. The plots show the absolute value of the differ- 
ence 8E N = E Q (N, D = const^-Eg 3 . The FSS is represented 
by a diagonal black line following the scaling from Eq. [8] All 
data sets initially follow this line. The FES regime corre- 
sponds to the various horizontal lines, where SEn saturates 
for different D to different values SEd that do not depend 
on N. In the inset we collect these values to show that they 
reproduce the expected behavior of the FES. The two regimes 
are separated by the appearance of a pronounced peak. Since 
we plot an absolute value, the peak is nothing more than the 
change of sign in the difference Eo(N, D) — E^ when moving 
from the FSS regime |8f to the FES regime |9f. b) The same 
plot for the HB model tells us that here the FSS is much more 
difficult to study, since all data-sets deviate very soon from 
the pure FSS prediction. 



E (D) -E (00) cx— . (9) 

where w = 2k and k is the same exponent of Eq. [2] and 
A is a positive non-universal constant. We show below 
that the same scaling holds also for MPS simulations of 
finite chains if N is big enough. This happens exactly in 
the FES regime defined by Q. 

The two regimes are very clearly distinguished in Fig. 
[2] where we present plots of the absolute value of the 
difference of the ground state energy obtained with MPS 
simulations and the exact value in the thermodynamic 
limit 



5E N , D = E {N,D)-E (oo) (10) 

as a function of N in a log-log scale. Note that in 
Fig. [2] we make an abuse of notation by using SEn — 
{5E NtD ) D=const . and SE D = {SE NtD ) N=const ., which can 
be only done as long as we specify what the constant 
value of D or N is. The data is collected from several 
simulations of the critical IS with PBC for chain lengths 
in the range 20 < N < 4 • 1 4 (panel a) and of the HB 
with PBC in the range 10 2 < N < 5 • I0 3 (panel b). D is 
going in both cases up to D = 64. Each line in the main 
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plot represents simulations performed for different N at 
fixed D. The FSS predictions of Eq. [8] are straight lines 
plotted in black. For small N, each set of data follows 
the prediction of Eq. |8j which is a clear signal of the FSS 
regime. The maximal N for which the FSS prediction 
holds increases with growing D as expected. However, 
each set deviates at some big enough N from the FSS 
prediction to eventually stabilize to a value of the energy 
difference that only depends on D. This is a clear foot- 
print of the FES regime scaling, as described in Eq. [9] 
In order to confirm this we have added to both panels 
insets where we have plotted several values of SEu for 
large fixed N as a function of D in a log-log scale. Simi- 
lar plots in the thermodynamic limit can be found in [29] . 
The linear fits (red lines) in the insets yield «js ~ 1.9776 
for N = 10 4 respectively k H b ~ 1.3025 for N = 3000. 
These values are compatible with the analytical values 
obtained for N -> oo in [10], namely n a ^ al ps 2.03425 
and Kfl l g l fa 1.34405 and thereby confirm the scaling of 
Eq.[9} 

Note that we are able to obtain much better precision 
for the IS than for the HB model at at the same com- 
putational cost . This is visible by comparing the panel 
a) to the panel b) and observing that for fixed D, the 
curves for the HB model deviate from the FSS at much 
lower values of N than the corresponding ones for the IS 
model. 



B. The transition between the two regimes 

In Figure [2] we can observe that for each line with con- 
stant D, the FSS region is separated from the FES re- 
gion by a well distinguishable peak in the absolute value 
of 6En. We would now like to show that this transition 
does not depend on the choice of the observable but that 
it indicates a global change in the wave function. 

To this end we can investigate the trace distance be- 
tween the exact ground state of a chain with N sites and 
the MPS obtained from a series of simulations with dif- 
ferent D. We have chosen the step size in D as small 
as possible, i.e. AD — 1. Since the exact ground state 
wave function is only available for very small systems 
due to the exponential scaling of the number of param- 
eters, we use as a reference state a MPS approximation 
of the ground state with very big D. For the N range in 
question, the biggest available bond dimension is D = 64. 
Note that the energy difference between the exact ground 
states and the reference states is much smaller than the 
difference to the MPS we want to compare to (see Fig. 

Figure [3] shows the trace distance between states with 
relatively small D and reference states for several differ- 
ent chain lengths N for both the IS and the HB models. 
Note that for every TV there is a jump in the trace dis- 
tance between states that are very far away from the 
reference state and states that are at least one order of 
magnitude closer to it. For the IS model the jump is very 




FIG. 3. (Color online), a) Quantum Ising model: trace 
distance between reference states with D\ — 64 and ground 
state MPS with bond dimensions D2 for several different chain 
lengths, b) Heisenberg model: trace distance between refer- 
ence states with D\ — 64 and ground state MPS with bond 
dimensions D2 for several different chain lengths. 

steep and for each line of constant N we can clearly iden- 
tify D[ as the biggest D in the left (FES) regime and D r 
as the smallest D in the right (FSS) regime. In this case 
the appearance of the jump evidently indicates the tran- 
sition from the FSS regime, where the trace distance is 
close to zero to the FES regime where the trace distance 
abruptly increases. For the HB model the transition is 
much smoother and we can not unambiguously define Di 
and D r for all lines with constant N. This happens pre- 
sumably due to the fact that taking MPS with D = 64 
as reference states is not accurate enough in the case of 
the HB model. 

However, at least for the IS model, we are now in the 
position to check if our intuitive expectation, that the 
transition occurs precisely when the correlation length of 
the finite size MPS reaches the size of chain as described 
in Eq. [5] is quantitatively correct. The correlation length 
of MPS with finite D reads according to Eq. [2] as £(-D) = 
k c ■ D K where k c is a proportionality constant. For our 
numerical study we obtain the parameters k c and k in 
the appendix [A] We can confirm that for each line of 
constant N the jump in the trace distance is consistent 
with our assumption, i.e. £(Di) < N < £,(D r ). 
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Furthermore we would like to mention that jumps also 
occur in other quantities at the same Di, for instance in 
the half-chain correlation function reported in [TH] (see 
figure 9 in that work for a plot of the jump for N — 500). 
The fact that the induced correlation length £(-D) in the 
FES is smaller than the size of the system, suggests that 
the state is completely unaware of the presence of the 
boundaries. This confirms our intuition that MPS in the 
FSS regime are faithful approximations of finite chains 
with PBC while MPS in the FES regime do not capture 
properties related to the boundary conditions. 

Summarizing, the main point of this section is that if 
one is interested in the effects of PBC, results collected 
for D smaller than D r should not be taken into account. 
Note that due to the residual dependence on D (see the 
appendix II B 1 for details on this point) one still has to 
extrapolate the results in the limit D — ¥ oo in order to 
obtain accurate results. If, on the other hand, one is 
interested only in local universal quantities (i.e. where 
boundary conditions are irrelevant), there is no point in 
simulating the system with PBC and one should rather 
perform a standard FES study [7|. 



1. The real scenario, a complex cross-over induced by 
corrections to the scaling 

We have seen above how in the extremal regions of Fig- 
ure^ the simulation results follow the behavior predicted 
by FSS respectively FES. In the intermediate region how- 
ever, the simulations display a behavior that can not be 
attributed to either regime. We would now like to point 
out that the real picture is somewhat more complex than 
the two-regime interpretation given above. 

The leading scaling behavior given in Equations [8] and 
[9] represents only the first terms of the series expansion 
of more complex analytic corrections. Thus these terms 
are accompanied by higher order terms called corrections 
to scaling. In order to understand the scenario we must 
consider the general Taylor expansion of a two variable 
function. Let us consider two variables A d and A n with 
the property that lim£)->.co A^ = and limAr_ i . 00 Ajv = 
0. Obviously these variables can be identified with the 
gaps proportional to the inverse of the correlation length 
induced by the system size N and by the finite matrix 
dimension D as defined in Eq. [T]and Eq. eq:feg. Part of 
the scaling ansatz consists in assuming that all universal 
quantities are universal functions of these two variables. 

Let us review the case of a one-scale problem. In this 
case, by neglecting higher than quadratic terms in the 
vanishing variable (e.g. A„) we get the following series 
expansion for some universal function g 



9a n = .9o + dAjvStoA 



1 



JV 



N 



(11) 



In the regime where A N <C 1, the first two terms are con- 
sidered the leading scaling behavior while the rest pro- 



vides only higher order corrections. If we now take a two 
scale problem 



,,An — fo,o + <9a d /ooA_d + 9ajv/ooAat 
+ l(di D fooA 2 D + di N f 00 A%) 



" ^AdAn/oqAdAat + • 



(12) 



in the regime where Ad <C A n <C A^ we are back to the 
previous situation and we can apply the one-scale ansatz 
of Eq. 11 to the function g(Ajv) = /(Ajv,0); the same 



thing is valid in the opposite regime A^r -C A D < A^ 
with the obvious substitution g{Ar)) = /(0, Ary). These 
two limits would correspond to what we have called in 
the main text the FSS regime and the FES regime (see 
Eq. [3]and[4]). 




FIG. 4. (Color online). Classification of MPS simulations 
of spin chains according to the simulation parameter pair 
{N,D}. All lines ending in the origin denote possible paths to 
approach the thermodynamic limit of critical systems when 
doing MPS simulations. 

Now in general, there is a huge regime where for ex- 



Aj, <C An < A 



In 



ample A 2 N <C Ad < Ajv or 
this case the leading scaling behavior is modified by cor- 
rections that are not proportional to the next power of 
the relevant variable but to the ratio among the two 
variables. Indeed if w e c onsider the scenario where 
< An in Eq. 12 we obtain 



A D < A N 



/ad.Ajv = /o,o+Ad <9a d /oo + dA N f{ 



'A 



D 



(13) 



How relevant the correction is clearly depends on the 
scale separation, i.e. on how close Ajv/Ad is to one. 
In the following we give a sketch of how this cross-over 
region looks like and we introduce two new terms: Finite 
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Entanglement- Size Scaling (FESS for) the region where 
the leading scaling is due to the finite size of the matrices 
and the corrections come from the size of the system 
and Finite Size- Entanglement Scaling (FSES) where the 
leading scaling is due to the size of the system and the 
corrections come from the size of the matrices. 

Figure [4] shows a classification of MPS simulations ac- 
cording to the simulation parameter pair {TV, D}. The 
thermodynamical limit can be approached by moving 
along any path towards the origin of the diagram {TV -1 = 
0, D^ 1 = 0}. However in order not to distort the scal- 
ing analysis by mixing the different TV and D related 
corrections, moving from one point to the next on the 
path should leave the ratio Ajy/Ajj unchanged. This is 
equivalent to the requirement that any path is completely 
determined by the path constant k = N/D K . 

We can distinguish three different regions and three im- 
portant lines in Fig. [4] In the region above the blue line 
which is defined by D = d N ^ 2 , the MPS bond dimension 
is large enough to represent the ground state exactly. Of 
course doing MPS simulations in this regime is point- 
less since the computational cost becomes exponential in 
TV and there is no advantage over exact diagonalization. 
Thus no matter which path towards TV — > oo we choose in 
this region, it is completely equivalent to FSS. The ma- 
genta line with TV -1 = represents the only path along 
which pure finite entanglement scaling (FES) holds. The 
red line represents the path along which the induced cor- 
relation length is equal to the system size, i.e. TV = £,(D). 
We will call this line in the following the critical line. For 
critical models without conformal invariance the critical 
line can be obtained using the method described in ap- 
pendix [5] Between this line and the FES line there is 
a region where TV > £(D), All simulations done in this 
region barely registrate the boundaries of the system and 
the fixed point MPS is more or less the same like that of 
a TV = oo simulation with same D. However there is a 
slight effect due to the finite size for points close to the 
TV = £(D) line as can be seen in Figure [l] This is why 
we call this region the finite entanglement-size scaling 
(FESS) region: the entanglement scaling predominates, 
but there is a small trace of finite size scaling behavior. 
The region between the critical line and the FSS-regime 
describes MPS simulations where £(D) < TV, which turn 
out to reproduce faithfully the long range correlations 
throughout the entire chain (see figure 9 in our previ- 
ous work 18J). The FSS aspect predominates in this 
region, however there is also the inherent error of MPS 
simulations with D < d N / 2 , so we call it the finite size- 
entanglement scaling (FSES) region. 

Despite the rigorous classification of regimes from Fig- 
ure [4] we will restrict ourselves in the following to dis- 
criminate merely between the regimes on different sides 
of the critical line. We do this in order to improve the 
readability. Thus we will refer to both FSS and FSES as 
FSS; analogously we will denote both FES and FESS as 
FES. 



N=1000 



a) 




b) 

FIG. 5. a) (Color online). Quantum Ising model: trace 
distance between ground state MPS with different bond di- 
mensions D\ and Di for a chain with TV = 1000 sites, b) 
Heisenberg model: trace distance between ground state MPS 
with different bond dimensions D\ and D2 for a chain with 
TV = 200 sites. By using relatively small bond dimensions we 
are able to localize the transition between the FES and FSS 
regime for each TV. This can be used in case we are interested 
in performing a FSS analisys by providing the lower bound 
D r such that the state we obtain is not orthogonal to the 
exact state. 



C. Minimal D for faithful simulations 

We can outline a direct application of the presence of 
a transition between a FSS regime and a FES regime. 
Suppose that we want to simulate a critical chain with 
PBC such that it is in the FSS regime, i.e. the properties 
due to the boundary conditions are faithfully reproduced 
since TV < £(D). In order to minimize the computational 
cost we would like to use the smallest possible D that 



captures the PBC topology. By looking at figure [3] it is 
clear that we would have to choose D = D r to this end. 
The problem is that in order to make that plot we had 
to use as reference states MPS with very large D = 64, 
which is exactly what we would like to avoid in this case. 
Fortunately it turns out that even without a large D 
simulation it is possible to detect the optimal D = D r . 
This is due to the fact that all MPS with D > D r have a 
much smaller trace distance among eachother than with 
MPS with D <D r . 

The trace distance among all states with D < 40 for a 
IS chain with N = 1000 is shown in figure [5ja). The plot 
is of course symmetric in D\ and D2 and we have omit- 
ted the points on the diagonal since they are trivially 0. 
The transition between the FSS and the FES regime is 
clearly distinguishable at the same location of the jump 
as in figure [3] but in this plot we used only MPS with 
relatively small bond dimension. Note furthermore that 
if D\ , D2 < D r the trace distance between these states 
is wildly oscillating. However if D\ and D2 are on differ- 
ent sides of the jump, profiles similar to figure [3] emerge. 
Now it is clear how we can find the optimal D = D r with 
the smallest computational cost possible: for a given N 
run the PBC simulations by increasing D in small steps, 
ideally AD = 1. After each simulation compute the over- 
lap with all previously obtained MPS and when the nice 
profile with the jump appears, we know we have reached 
D = D r . The same strategy can be employed for the 
HB model, however, just like in figure [3] the transition is 
much smoother in this case. 

As a side remark note that due to the fast decay of 
the eigenvalues of the MPS transfer matrix one can com- 
pute the overlaps with computational cost scaling like 
0(nD 3 ). The meaning of n and the method how to 
achieve this is described in 1181. 



D. Thermodynamic limit of the transition 

What can figure [3] tell us about the behavior of the 
transition between the FSS and FES in the limit N — > 

00? 

For the IS model, qualitatively the height of the jump 
seems to remain constant for increasing N. The trace 
distance between MPS with bond dimensions Di and D r 
to the reference state also seems to remain more or less 
stable but this is of course not enough evidence for the 
persistence of the transition in the thermodynamic limit. 
In appendix [B] we present a detailed analysis that shows 
that for the IS model i) the N — > 00 limit of the trace dis- 
tance between the exact ground state (approximated by 
a reference state) and MPS obtained in the FSS regime 
is strictly bigger than zero, ii) the same limit for the 
trace distance with respect to MPS obtained in the FES 
regime is zero, ii) implies that states in the FES regime 
are globally orthogonal to the exact ground state of the 
PBC chain. As we already mentioned above this does 
not affect the possibility to extract local universal infor- 



mation from those states. However ii) clearly shows that 
MPS in the FES regime are globally not a good approx- 
imation for the ground state of the IS model with PBC. 

Unfortunately we cannot obtain the same conclusions 
for the HB model. Presumably this is due to the fact that 
the reference states that we use are not a good enough 
approximation of the true ground state of the model in 
this case. This becomes clear if we look again at figure [T] 
for the IS model the D = 64 states have a much better 
precision than the MPS we compared them to in order 
to prove the persistence of the transition in the thermo- 
dynamic limit (see Appendix [B] for details). For the HB 
model on the other hand, the D = 64 line covers almost 
three orders of magnitude in the relative precision plot; 
at its maximum it is over one order of magnitude above 
the points belonging to MPS that we must compare the 
reference states to in order to perform our analysis of the 
thermodynamic limit (e.g. the data points with N = 100 
and D = 48). The D = 128 line in right plot of figure [l] 
seems to fulfill similar requirements like the D = 64 line 
in the left plot. However, in that regime, for N -C £(D), 
the PBC algorithm is very inefficient and it would take 
unreasonably long to obtain the data points for D = 128. 



E. The scaling function 

Finally we conduct an analysis of the scaling of MPS 
simulations across the entire interval N/£(D) € (0, 00) 
which covers all possible pairs {N, D}. This is very much 
in the spirit of the scaling analysis performed by Nishino 
et al. for classical 2D systems in Ref. [27]. The main dif- 
ferences are that in our case the energy difference SEn^d 
can take both positive and negative values, and that we 
obtain the effective correlation length £(D) from an anal- 
ysis of the humps in the relative precision of the energy 
(see Appendix [A] for details) instead of using the ratio 
between the two biggest eigenvalues of the MPS transfer 
matrix. 

Analogously to Nishino, we first eliminate the FSS scal- 
ing from |5Sjv,d| and then we plot the result (in our case 
this is N 2 \SE N ^ D \) as a function of N/£(D). The fact 
that all data (with exception of the D = 64 points for 
the IS model) collapses into a single curve justifies the 
assumption that 



SE N>D = E (N, D) - £ (oo) = f{N/ ^ 2 {D)) • (14) 

with some scaling function f(x) that is not exactly 
known. What we can easily write down however is its 
asymptotic behavior 



lim /(x) = 

x— yQ 

lim f(x) = 
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FIG. 6. (Color online), a) Quantum Ising model: log-log 
plot of iV 2 |<5£^ d\ versus N/£(D) that illustrates the collapse 
of the data into a single curve. The points with D — 64 
slightly deviate from the curve traced by data points with 
smaller D. b) Heisenberg model: log-log plot of N 2 \SEn,d\ 
versus N/£(D) that illustrates the collapse of the data into a 
single curve. 



For the IS model we have used for ((D) the expres- 
sion obtained from the hyperbola fit in Fig. [8] of the Ap- 
pendix |Aj i.e. ((D) = 3.810 • D 2 042 . Note that in plot 
a) of the Fig. |6j the data for different D collapses al- 
most perfectly in the extremal regimes N <§C ((D) and 
N > ((D). There is a slight deviation of the D = 64 
curve that can be explained if we look at figure 1 in [21] 
(there the D = 64 data point also slightly deviates from 
the line that is traced by the points with D < 64). In 
the regime where N s» ((D) the curves do not collapse so 
nicely which is a manifestation of the fact that the humps 
in Fig. [I] are so different for the IS model. 

For the HB model we have used for the effective cor- 
relation length ((D) = 3.647 • D 133S as obtained from 
the hyperbola fit in Fig. |8j Plot b) of Fig. [6] shows an 
almost perfect collapse even in the region N s» ((D). 
Presumably this is due to the fact that for the HB model 
the humps in Fig. [T] are much more similar among each 
other than in the case of the IS model. 



III. CONCLUSIONS 

An accurate analysis of MPS simulations of critical 
spin chains with PBC reveals the appearance of two 
regimes. The FSS regime where the energy gap of the 
system is induced by the size of the system and the FES 
regime where an effective energy gap is induced by finite 
D. While in both regimes local universal quantities can 
be extracted by studying the scaling of the observable 
with respect to the relevant variable (the size of the sys- 
tem for the FSS or the size of the MPS matrices for the 
FES regime), we have shown that for the Quantum Ising 
model, states in the FES regime are orthogonal to the ex- 
act ground state in the thermodynamic limit. Intuitively 
this happens due to the fact that for MPS simulations 
in the in the FES regime, the induced correlation length 
is smaller than the system size and thus the MPS is not 
aware of the size of the system. Since in critical systems 
the boundary conditions strongly affect global properties 
of the system, this result seems quite natural. 

Our results can be interpreted as a further bench- 
mark for recently introduced algorithms that try to lower 
the computational cost of PBC simulations with MPS 
[T51 (see Appendix [C| . Here we provide strong 

hints that in order to correctly describe the ground state 
of a finite chain with PBC for critical systems, these al- 
gorithms should be used with care in order not to obtain 
wave functions that are orthogonal to the exact ones. 
What one would indeed interpret as the MPS tensor for 
a PBC chain, in some regime could turn out to be closer 
to the MPS tensor of an infinite OBC system. 

However, considering that for OBC systems the ap- 
proach to the thermodynamic limit is by no means slower 
than for PBC systems (for both the ground state energy 
converges to the thermodynamic limit as a function of 
the appropriate correlation length like (~ 2 ), our results 
can be also used in a constructive way. In order to ex- 
tract universal information about local operators, one is 
better off by using FES rather than FSS, since simula- 
tions in the FES regime have a much better scaling of 
the computational cost. 

Things are more complex if one is interested in global 
observables, such as e.g. two point correlation functions 
at half chain length. For PBC systems the scaling anal- 
ysis must be performed in this case on paths with con- 
stant k = N/D K that lie completely in the FSS regime. 
The computationally least expensive such path is the one 
where for every given N, the MPS bond dimension D is 
just big enough such that (jj > £jv- We have shown 
how that minimal D can be found for any N by looking 
at the overlap between MPS with increasing D until the 
discrete transition between the FES and the FSS regime 
is detected. Regarding the scaling exponent K, we have 
been able to numerically confirm the theoretically pre- 
dicted values with an accuracy of approximately 0.4% 
for both the Quantum Ising and the Heisenberg models. 
Furthermore we have shown in Appendix |A 1| how the 
analytical expression for k, originally derived in |10j . can 
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be obtained in an alternative way. 

Following Nishino's analysis for 2D classical sys- 
tems [57] we have shown that also for MPS simulations 
of ID quantum systems the scaling of the MPS ground 
state energy in simulations with finite iV and D obeys 
a two-parameter scaling function. Finding an analytical 
expression for this function is something that still has to 
be done. 

A further interesting future line of research is to un- 
derstand how to extract information about the operator 
content of the Conformal Field Theories related to the in- 
frared behavior of the studied critical spin systems (that 
strongly depend on boundary conditions 0[26]) directly 



out of the MPS tensors. 
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been observed that the usage of finite D in MPS simula- 
tions close to the critical point leads to an effective shift 
of the critical point (see Fig. 2 in |? ). This observation 
led us to the idea of equating the corrections in Eqs 
and ( A2 1 and identifying £, 
= k, ' 



assumption £j 



with £rj. 
D K this yields 



(Al) 



Together with the 



P r {b,D) 



A! ■ D~ 



(A3) 



where we have collected all constants into A' — A/{k c -0). 

In the large D limit (required due to our assumption 
of working in the scaling limit), the residual probability 
reads according to [TU] 



P r (b,D) 



2be~ b D 
\ogD- 2b ( 



-(log D) 2 /4b 



(A4) 



where 



Appendix A: Effective correlation length 

1. Analytical results 

Recently it has been shown numerically that any MPS 
simulation of an infinite spin chain leads to the emergence 
of an effective correlation length induced by the finite 
rank D of the MPS matrices, even if the studied system 
is critical [7J . In Ref. [TU] the authors relate the numerical 
observation that £(~D) ex D K from Ref. [7] to analytical 
results on the spectrum of the MPS transfer matrix [31] 
and to well-known results from conformal field theory [5J 
|2"8] in order to derive an analytical expression for the 
exponent k. Here we derive the same results in a different 
way. 

The starting point for our argument is the same like 
the one in [10] , namely that corrections to the exact 
ground state energy in the thermodynamic limit can 
have different origins. On one hand conformal invari- 
ance yields in the vicinity of the critical point (i.e. e = 
|A - \crit\/Krit < 1) according to Refs. [51 1101125] 



£o(&) = £0(00) + 



.4 

if 



(Al) 



where A is a non-universal constant. On the other hand, 
MPS simulations with finite D yield according to Refs. 



Eo(Zd) =£0(00) 



P r (b,D) 



(A2) 



where (3 is a non-universal constant, P r (b, D) is the resid- 
ual probability due to the usage of finite D and b is re- 
lated to the dominant eigenvalue of the reduced density 
matrix of the half-chain (see Refs. [HUGH]). Now it has 



(A5) 



and c denotes the central charge in the associated con- 
formal field theory. Inserting (A4) and (A5) into (A3) 
yields after several steps 



CH 



6 — CK 



A' ■ D~ K 



(A6) 



Equating the exponents in ( A6 ) yields a quadratic equa- 
tion for k with the solutions 



k± = 



(A7) 



The physical root is the one that is positive for all values 
of c, i.e. 



(A8) 



which is exactly the result obtained in Ref. [TU] . 



2. Numerical results 

In this appendix we show how the effective correla- 
tion length £(D) emerges in our simulations of finite spin 
chains with PBC. As the scaling of the algorithm [T5] is 
quasi-independent of the chain length N we can use it to 
approximate ground states of arbitrary long chains with 
PBC. The relative precision of the MPS ground state en- 
ergy for a given D is plotted in figure [T] as a function 
of N. Each of the lines contains a hump which can be 
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FIG. 7. (Color online). Quantum Ising model data. Left: relative precision of the PBC-MPS ground state energy (blue/dark) 
and relative precision of the state used as an input for the PBC algorithm (i.e. obtained by inserting the TL-MPS, cyan/bright) 
as compared to the exact result for D = 8 and D = 16 (inset). Right: ratio between the PBC-MPS energy and the one of the 
input state (TL-MPS) for D = 8 and D — 16 (inset). Fitting a degenerate hyperbola in form of two straight lines yields a well 
defined point (the intersection point) whose value as a function of D is proportional to the effective correlation length. 



interpreted as the evidence for a finite correlation length 
£(£>). In order to see this let us have a look at how the 
hump emerges. The left part of figure [7] shows a com- 
parison between the relative precision of the PBC-MPS 
ground state energy (i.e. the MPS towards which the al- 
gorithm in |18j has converged) and the relative precision 
of the energy for the MPS that we had used as a starting 
point for the gradient search. As explained in [18] this 
is the local MPS tensor obtained by imaginary time evo- 
lution 29J for a chain in the thermodynamic limit (TL) 
when it is used in the finite PBC geometry. One can see 
that for a given D, on the left side of the hump there 
is considerable improvement in the precision of the en- 
ergy between starting and ending point of the gradient 
search. As one approaches the hump from the left, the 
improvement decreases in order to vanish completely on 
the right side. This can be interpreted as follows: if N 
is too large for a given D, the finite chain looks for a 
local MPS-tensor as if it would be infinite. Sites that lie 
further apart than a certain correlation length £(D) effec- 
tively do not see each other. The transition to this region 
happens more or less smoothly since for growing powers 
of the MPS transfer matrix T, the subspace spanned by 
these powers gets smoothly restricted to the dominant 
eigenvector i.e. T N \ N:g> ^ D) « Af |Ai) (A x |. 

Thus the humps must represent some evidence for the 
emergence of a finite correlation length, but how can we 
extract some reliable numbers from them, as they differ 
considerably in shape and width? The answer is given by 
the right part of figure [7J We have observed empirically 
that if we make a log-log plot of E f mal (N) / E l mtial (N) 
[32) we obtain approximately two straight lines connected 
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FIG. 8. (Color online). Quantum Ising (upper) and Heisen- 
berg model (lower): linear fit of the logarithm of the effective 
correlation length as a function of the bond dimension D. 

by a small piece that is more or less smooth. This picture 
is reminiscent of a rotated hyperbola. We know further- 
more that in the large ./V limit all points have ordinate 
0. This suggests to fit a hyperbola that is degenerated to 
two straight lines through our data. The intersection of 
these lines is a well defined point which should be pro- 
portional to £,(D). 

Figure [8] shows log-log plots of the effective correla- 
tion length as defined above for both the Quantum Ising 
and the Heisenberg models. After fitting straight lines 
through each of the data sets we can read off the scal- 
ing £(£>) = k c -D K with {k w 2.042, k c w 3.810} /s and 
{k w 1.338, k c « 3.647}#s. Comparison with the ana- 
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lytical results (i.e. K a ^ al « 2.03425 and K a ££ l « 1.34405) 
yields a difference of roughly 0.4% for the Quantum Ising 
model and of roughly 0.43% for the Heisenberg model. 
These results are the ones we refer to in Sec. Ill Bl as the 
ones fulfilling £(£>;) < N < (,{D r ). 

An alternative way to extract the effective correlation 
length is obtained by interpreting the abscissa of the min- 
imum of each curve in Fig. [2] as a length proportional to 
£(D). Fitting a straight line through these minima in 
a log-log plot of N(D) yields for the IS model the ex- 
ponent kis ~ 2.0293 which approximates the analytical 
result with an accuracy of roughly 0.24%. On one hand 
this result is closer to the analytical value than the one 
obtained using the degenerated hyperbola fit. On the 
other hand, if we want to predict the bond dimension D 
for which the jump in the trace distance occurs in simu- 
lations with fixed TV, it turns out that the value obtained 
using this fit in does not always coincide with the actual 
values observed in figure [3] As mentioned above, the 
degenerated hyperbola fit satisfies this consistency test, 
which is why we prefer using that method to extract an 
approximation for k. Furthermore the plots in figure [2] 
require knowledge of the exact ground state energy in the 
thermodynamic limit, which is not always available. The 
strategy with the hyperbola fit on the other hand does 
not require any analytical results and thus can always be 
used. 



Appendix B: Detailed treatment of the 
thermodynamic limit of the transition 

In this appendix we present the details we used for the 
conclusion drawn in section III Dl of the main text. As 
mentioned above in appendix |IIB 1| a reliable analysis 
of the thermodynamic limit can only be made properly 
if we move towards it on paths of constant k = N/D K . 
However this analysis provides conclusive results only for 
the IS model which is why we skip presenting the results 
obtained for the HB model. As mentioned in the main 
text, the reason why this method fails for the HB model 
is that the reference states are in that case not precise 
enough. 

As a first step let us normalize the tensors in our states 
such that the largest eigenvalue of the MPS transfer ma- 
trix T is equal to one (i.e. Ai = 1 and Xi > Xj, Vz < j). 
This yields for the norm of such a state 



(*(£>, N)\*(D, AO) = Tr(T ) = 1 + V Af (D, N) 



■£■ 

=2 



(Bl) 

We will always use in the following lower-case greek let- 
ters to denote states that are normalized to one and 
upper-case letters for the corresponding state normalized 
according to (Bl), i.e. 



I*) 



(B2) 



For the computation of the trace distance between refer- 
ence states and states lying on a curve with fixed k we 
need the absolute square of the overlap which becomes 



mD k , N ,N)\^(D iei ,N))f 



\(V(D k , N ,N)\*(D I<iU A0)| 5 



(*(D ktN , N)\V(D k>N , AT)) (*(Arf, N)\H>(D xei , N)) 







2 


1 + ES W AfpMV,A0 







(B3) 

In the numerator we have used Hi(D k ^, D TC f, N) =: 
fJ>i(k, N) to denote the eigenvalues of the overlap transfer 
matrix 



T ovlp (k,N) 



d 

£ 

i=l 



Ai(D k . N , N) ® A*(D ie f, N) , (B4) 



where the Ai(D, N) represent as usually the matrices of 
a translationally invariant MPS with N sites and vir- 
tual bond dimension D. Similarly we will use for the 
eigenvalues of the MPS transfer matrix the notation 
Xi(k,N) := Xi(D k ,N , N) in the following. This can be 
done since we need only two of the quantities (D,N,k) 
to uniquely specify the point of the phase diagram that 
we want to refer to. 

The crucial argument in favor of the persistence of 
a discrete transition between the two regimes in the 
thermodynamic limit will be the fact that in this limit 
Hi{k,N) converges much faster to 1 in the FSS regime 
(i.e. for k < k c ) than it does in the FES regime (i.e. 
for k > k c ). In fact we will show below that in the first 
case limAr^oo (k, N) = 1 while in the second case we 
have lmijv->co Mi i^j ^f) = 0- The other contributions in 
the numerator of (|B3|) will turn out to be negligible for 



N — > oo, i.e. limjv_ s . 00 {k, N) = for any k and all 
i > 1. Furthermore we will show that the denominator 



of (B3) remains finite in all cases such that we will be 
able to conclude that the overlap of the quasi-exact |33j 
ground state with states in the FES regime converges to 
zero in the thermodynamic limit. Along the same lines 
we will argue that the overlap of the quasi-exact ground 
state with states in the FSES regime is always larger than 
zero in the thermodynamic limit, thereby concluding that 
a detectable transition between the two regimes persists 
for N — > oo. 

To this end we have considered three paths in the FSS 
regime (k sa 0.37,0.54,0.97) and two paths in the FES 
regime (k ss 18.0, 58.7). The exact data points (N, D) for 
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four of these paths are listed in table [TJ Note that since 
TV, D E N the exact value for k — N/D K varies slightly 
within each path. 



2.03425 



k w 0.37 < k c 


k « 0.54 < fc c 


k f 


a 18.0 > k c 


k r 


a 58.7 > fe c 


k 


TV 


D 


k 


TV 


D 


k 


TV 


D 


k 


TV 


D 


0.374 


122 


17 


0.538 


118 


14 


17 


9 


300 


4 


58 


9 


1000 


4 


0.373 


206 


22 


0.540 


198 


18 


17 


8 


470 


5 


59 





1580 


5 


0.371 


288 


26 


0.540 


298 


22 


18 


3 


700 


6 


58 


7 


2280 


6 


0.371 


386 


30 


0.536 


384 


25 


18 


1 


950 


7 


58 


8 


3130 


7 


0.369 


526 


35 


0.539 


560 


30 


18 


2 


1250 


8 


58 


6 


4100 


8 


0.369 


690 


40 


0.537 


810 


36 


17 


7 


1550 


9 


58 


6 


5210 


9 


0.368 


1000 


48 


0.535 


1000 


40 


18 





1950 


10 


58 


5 


6450 


10 














17 


9 


2350 


11 


58 


4 


7830 


11 














17 


9 


2800 


12 


58 


4 


9350 


12 



TABLE I. Data points constituting several of the investigated 
paths with roughly constant k depicted in figure [TT] 




FIG. 9. (Color online). Quantum Ising model: scaling of 
the eigenvalues of the overlap transfer matrix for five paths 
with roughly constant k — N/D K in different regimes. The 
exact pairs (TV, D) for the data points are given in table [I] In 
the legend we have only explained in detail what the differ- 
ent markers mean for the path with k S3 0.37. The markers 
for the other paths follow the same pattern. The full red 
lines are linear fits through the data for each f-ii(k, TV) respec- 
tively. The legend entries for these lines contain the values 

H?i(*),ai(*))- 

Let us first investigate how the numerator of the ratio 
( |B3[ ) behaves. We have observed that if we look at the 
eigenvalues /ii(k,N) along paths with constant k, then 
1 — fj,i(k,N) scales polynomially in TV as can be seen in 
the log- log plot of figure [9j such that we have 



IH(k,N) = 1 



TVft« 



(B5) 



Figure [9] shows a log- log plot of 1 — TV) for all k and 
fixed D re f = 64. The numerical values of aj(fc) and ft(fe) 



k = 2.03425 





k « 0.37 < fe c 


k S3 0.54 < k c 


k S3 0.97 < k c 


i 


ft 


OLi 


Pi 


a% 


Pi 




1 


4.06477 


57128.2 


3 


49773 


3170.6 


2.90998 


187.3 


2 


0.66454 


0.14622 





64177 


0.12814 


0.60463 


0.10416 


3 


0.51554 


0.19106 





51270 


0.19112 


0.48326 


0.16400 


4 


0.58048 


0.37486 





57835 


0.37607 


0.55326 


0.33117 


5 


0.48173 


0.26650 





51001 


0.31294 


0.51875 


0.32733 


6 


0.50660 


0.35546 





50828 


0.36216 


0.49502 


0.33813 


7 


0.47451 


0.31883 





46129 


0.30227 


0.42975 


0.26069 


8 


0.46673 


0.35474 





46928 


0.36415 


0.45436 


0.33884 


9 


0.48042 


0.42354 





47013 


0.41012 


0.44586 


0.37264 


10 


0.51091 


0.55418 





50601 


0.54819 


0.48431 


0.49597 



TABLE II. Scaling of fii(k,N): parameters Pi(k) and ai(k) 
for the fitting of 1 - /^(fe, TV) = a i (fc)/TV /3i(fe) for paths in the 
FSES regime. 



K = 2.03425 





k S3 18.0 > fee 


k S3 58.7 > fe c 


i 


ft 




ft 


OLi 


1 


0.94079 





20604 


0.95736 


0.46203 


2 


0.72904 





74874 


0.75036 


1.14015 


3 


0.46762 





25046 


0.43380 


0.23299 


4 


0.42154 





26209 


0.49195 


0.62481 


5 


0.37208 





23363 


0.37851 


0.33497 


6 


0.38713 





38566 


0.44232 


0.87896 


7 


0.36770 





38856 


0.38792 


0.65807 


8 


0.40633 





58072 


0.41378 


0.94025 


9 


0.42088 





80086 


0.45836 


1.67698 


10 


0.40215 





80529 


0.43306 


1.58778 



TABLE III. Scaling of /i;(fc,TV): parameters ft(fc) and on(k) 
for the fitting of 1 - /Li* (ft, TV) = a i (fe)/TV /3i(fc) for paths in the 
FESS regime. 



for the 10 largest Hi(k,N) are listed for the paths in the 
FES regime in table [IT] The equivalent data for paths in 
the FES regime can be found in table |III| We see that 
in the FSS regime for i = 1 we have f5\{k) > 1 while in 
all other cases we get 0i(k) < 1. This means that for 
TV — > co the overlap ( B3 1 always converges to zero in the 
FES regime due to 



lim (1 

N— s-oo 



TV^ 



) iV =0 V/3<l,a>0 



(B6) 



and due to the fact that the denominator is always larger 
than zero (in fact it is always larger than one). In the 
FSS regime on the other hand, the i = 1 terms in the 
numerator of ( B3 ) survive in the thermodynamic limit 
due to 



N 



lim (1 



a . 
TV?- 



N 



1 V/3>l,a>0 



(B7) 
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K = 2.03425 




k = 2.03425 




k « 0.37 < fc c 


k « 0.54 < k c 


k w 0.97 < k c 




k w 18.0 > k c 


k « 58.7 > k c 


i 






ft 


Oti 


ft 




i 


ft 


on 


ft 




2 


0.60788 





09304 


0.65085 


0.12184 


0.62045 





09914 




2 


0.84555 


0.42728 


0.79472 





28481 


3 


0.41227 





08847 


0.46692 


0.12433 


0.44144 





10401 




3 


0.76096 


0.82750 


0.67333 





41376 


4 


0.43118 





12846 


0.48404 


0.17837 


0.45909 





14987 




4 


0.76823 


1.11696 


0.68922 





59758 


5 


0.30225 





07647 


0.36915 


0.11385 


0.35609 





10226 




5 


0.61161 


0.54136 


0.51839 





26161 


6 


0.34175 





10783 


0.40765 


0.15889 


0.39654 





14462 




6 


0.55554 


0.41101 


0.49502 





25624 


7 


0.40046 





17074 


0.44681 


0.22606 


0.43958 





21126 




7 


0.42416 


0.19075 


0.52876 





42488 


8 


0.41688 





19910 


0.46402 


0.26460 


0.45243 





24089 




8 


0.36840 


0.13854 


0.42626 





21306 


9 


0.33412 





13003 


0.36584 


0.15862 


0.36051 





15116 




9 


0.42172 


0.22169 


0.40383 





19455 


10 


0.28957 





10620 


0.34067 


0.14540 


0.33743 





14045 




10 


0.39138 


0.19807 


0.31502 





10918 



TABLE IV. Scaling of \i(D xei ,k, N): parameters ft(fc) and TABLE V. Scaling of \i(D Icf ,k,N): parameters ft(fc) and 
Oi(ft) for the fitting of 1 - \i(D Tel , k, N) = a i {k)/N l3iW for on(k) for the fitting of 1 - Xi(D rel ,k,N) = a i (k)/N l3i ^ for 
paths in the FSES regime. paths in the FESS regime. 



However this is not enough in order to show that the 
overlap is strictly larger than zero in this regime. A di- 
verging denominator in the limit N — > oo could spoil this 
line of reasoning, so we have to convince ourselves that 
both factors in the denominator of ( B3 ) remain finite in 
the thermodynamic limit. 
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FIG. 10. (Color online). Quantum Ising model: scaling 
of the eigenvalues of the MPS transfer matrix for D rc { = 64 
and all N occuring in table [I] In the legend we have only 
explained in detail what the different markers mean for the 
N in the path with k m 0.37. The markers for the other paths 
follow the same pattern. The dotted red lines are linear fits 
through the data for each \2(k,N) respectively. The legend 
entries for these lines contain the values (— ft>(fc), oi2(k)). 

Let us first treat the norm of the reference MPS since 



this turns out to be the easier one. Figure 10 shows a log- 
log plot of 1 - Xi (k, N) for all k and fixed D ref = 64. The 
numerical values for i < 10 are given in tables [TV| and |V| 
For large chains with N > 1000 figure [TO] clearly indicates 
polynomial scaling in N. Note that for small chains with 
N < 1000 the plot deviates from the nice linear behavior 
that we see for TV > 1000. The reason for this are numer- 
ical errors in the computation of the ground state MPS. 



This effect can also be seen in figure [T] for N < 1000 the 
algorithm we use cannot minimze the energy beyond a 
relative precision of roughly 8 • 10 -11 even if we decrease 
N while keeping a constant D — 64. Apart from that, 
the fitting in figure 10 yields all ft(fc) < 1 for i > 2 thus 
we can conclude that the norm (ty(D re f, iV)|^E'(Z3 re j, N)) 
converges to one in the thermodynamic limit when we use 
the normalization prescription (Bll. 




FIG. 11. (Color online). Quantum Ising model: scaling 
of the eigenvalues of the MPS transfer matrix for five paths 
with roughly constant k — N/D K in different regimes. The 
exact pairs (N, D) for the data points are given in table [T] In 
the legend we have only explained in detail what the differ- 
ent markers mean for the path with k w 0.37. The markers 
for the other paths follow the same pattern. The full red 
lines are linear fits through the data for each \?(k, N) respec- 
tively. The legend entries for these lines contain the values 

(-&(*), 02(fc)). 

The norm of the states \^(Dk,N> N)) along paths with 
constant k also turns out to converge to a finite value even 
though the argument is a bit trickier in this case. The 
scaling of the largest eigenvalues \(k,N) for each path 



is shown in figure 11 The numerical values for i < 10 are 
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k = 2.03425 





k w 0.37 < fe c 


k « 0.54 < fc c 


k « 0.97 < fe c 
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ft 
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1.37029 
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5.00085 
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5.42433 
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0.99397 


6.19275 





98129 


6.13836 





97395 


6.67115 
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1.00992 


11.36417 





97979 


10.64645 





96439 


11.79694 


6 


0.99749 


12.14274 





96755 


11.35451 





95325 


12.79517 


7 


0.99290 


15.20975 





96478 


14.29878 





94933 


15.72285 


8 


0.99385 


16.84721 





94982 


14.52223 





93257 


16.11736 


9 


1.03853 


23.13857 





98133 


19.63216 





94593 


20.59742 


10 


1.03720 


25.64027 





95765 


19.24667 





92658 


21.15268 



This allows us to approximate the overlap (with a quasi- 
exact state) towards which MPS simulations in different 
regimes converge to (on the paths we considered) as 



JV 



lim (ip(k,N)\iP(D iei ,N)) k<k «0.45 



lim (^(k,NM(D Tcf ,N)) k>kc =0 



(B10) 



Thus we can conclude that the thermodynamic limit of 
the overlap in the FSS regime is always greater than zero 
proving that there is indeed an discrete transition from 
the FSS regime to the FES regime where the overlap 
becomes zero. 



TABLE VI. Scaling of Xi(k,N): parameters ft(fc) and oa(k) 
for the fitting of 1 - Xi(k, N) = « i (fc)/Af 3i(fc) for paths in the 
FSES regime. 
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k w 58.7 > k c 
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ft 
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ft 
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1.01831 


40.25148 
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01731 


135.41657 


3 


0.98774 


81.55159 





98683 


264.69533 
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0.97128 


88.78991 





97053 


282.91558 


5 


0.92704 


111.61975 





92647 


337.75755 


6 


0.88126 


94.42375 





88080 


270.69784 


7 


0.88252 


104.24989 





88232 


299.97710 


8 


0.77072 


64.57797 





77072 


162.76623 


9 


0.79065 


82.78179 





78996 


212.48337 


10 


0.76184 


72.53475 





76185 


180.90762 



TABLE VII. Scaling of Ai(k,N): parameters ft(fc) and cti(k) 
for the fitting of 1 - Aiffc, N) = a i (fc)/7V' 3i(fc) for paths in the 
FESS regime. 



given in tables VI and VII We see that most of the ft (k) 
are very close to one for small i in contrast to the values 
obtained for Xi {D re f , N) which are all well below one. In 
fact some of the ft(fc) are even bigger than one sug gesti ng 

^ of 



Bl 



that liniAr_ i . 00 Af^ = 1 in these cases. In section 
this appendix we give evidence for the fact that even 
if ft(fc) > 1 in some cases, the number of these values 
remains finite for any k. Furthermore we argue that in 
these cases it is reasonable to assume that we actually 
have ft(fc) = 1 which yields in the thermodynamic limit 



lim Xf(k,N) 

iV— »-oo 



lim (l--i) w = exp(-a i ) . (B8) 

iv— >oo iv 



Summing up all relevant contributions then yields for the 
norm of states in the different regimes 



JV 



lim (^(k,N)\t>{k,N)) k 



k<k c 



2.2 



hm {Hf(k,N)\9(k t N))k>k c < 2 - 



(B9) 



1. Scaling of Xi(k,N) 

The first ten parameters at and ft for the MPS transfer 
matrix eigenvalues Xi(k, N) on paths in the FSES regime 
are given m ta ble [yi] the ones for paths in the FESS 
regime m table |VII| In the FESS regime we have ft > 1 
which then yields a contribution of limjv_>oo ^2 (^) = 1 
to the norm. For i > 2 we clearly see how the ft rapidly 
decay below one, thereby making sure that the corre- 
sponding contributions to the norm become zero in the 
thermodynamic limit. This means that if we approach 
the thermodynamic limit on paths in the FESS regime 
and always normalize the MPS according to (Bl I, i.e. 
Ai = 1, the norm of these states does not get bigger 
than two. In fact it is very likely that the true contri- 
bution of A2 is de facto zero [33]: for TV as big as 10 9 , 
using the values for ct2 and ft given in table |VII[ we get 
A^(fc = 18.0) sa 2 • 10~ 12 and Af (k = 58.7) « 4 • lO" 42 . 

In the FSES regime on the other hand the ft seem to 
oscillate randomly around one so we must look at the be- 
havior of larger i in order to see if and when they decay 
below one, which is what we ultimately need in order to 
show that the norm of these states remains finite in the 
thermodynamic limit when the normalization prescrip- 
tion ( Bl I is employed. 

Figure [12] shows a log-plot of the first 200 ft in the 
FSES regime and of the first 120 in the FESS regime. All 
curves are approximately straight lines in this plot which 
means that the ft decay exponentially with i. Remember 
that on the paths we chose to investigate in the FESS 
regime, the MPS with the largest virtual bond dimension 
have D = 12, thus we cannot fit any parameters ft for 
i > 121 since we have only one data point available there. 
For 100 < i < 121 we have only two data points, namely 
the ones for D = 11 and D = 12 (see table [i] in the 
main text) which is usually not the best premise for an 
accurate fit. Nonetheless the ft fitted in this range obey 
the same exponential decay observed for smaller i, where 
more data points are available. The inset in figure [i~2] 
shows a zoom into the region with i < 20. While for i < 8 
all ft in the FSES regime are very close to one, we observe 
that for larger i, the k = 0.37 line is visibly above the ft = 
1 line. This would suggest that in the thermodynamic 
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200 



FIG. 12. (Color online). Quantum Ising model: log-plot of 
the first values of @i(k) for i < 200. A zoomed view on the 
first 20 values is shown in the inset. 




the norm in the thermodynamic limit solely depends on 
on due to 



lim A 



N- 



(k, N) = 1 lim 

N-hx. 



(l-^) N = eX p(-a i ) . 

(BH) 

Figure 13 shows the behavior of exp(— a*j) for the paths in 
the FSES regime and % < 60, which according to figure 12 
is the problematic z-range. We see how all contributions 
rapidly decay below machine precision. Note that the 
black line (i.e. the path with k = 0.97) is for i > 17 sev- 
eral orders of magnitude above the k = 0.37 and k = 0.54 
lines, which is due to the fact that the corresponding Pi 
are so much smaller than one in this region, that the as- 
sumption f3i « 1 simply does not hold, and the actual 
contribution to the norm converges to zero. Note fur- 
thermore how for small i all three lines are almost on 
top of each other meaning that the values to which the 
norm converges in the thermodynamic limit for MPS on 
different paths in the FSES regime will be very similar. 
In fact we get 



lim N)\V(k, A0) fc =o.37 = 2.261646939734277 



JV 

lim (*(fc,AOI*(fc,AO) fc=0 .54 

iv— >-oo 



2.236037631274709 



lim iV)|*0, N)) k=Q g7 = 2.225635928039641 

iV— »oo 

(B12) 

which completes our argument that the norm of the MPS 
remains finite on any path in the FSES regime. 



FIG. 13. (Color online). Quantum Ising model: contribution 
of the eigenvalues with /3i(k) w 1 to the norm of the states in 
the thermodynamic limit, i.e. limjv-).oo (fe) « exp(— oti(k)). 
Note how all contributions fall off exponentially below ma- 
chine precision at i « 14. 



limit the eigenvalues Ai>8 would each yield a contribution 
equal to one to the norm while the contribution from the 
Aj with i € {2,3,4,6,7,8} would vanish, since in these 
cases Pi < 1. This makes however no sense since the A^ 
are decreasingly ordered, i.e. Aj > Xj if i < j. This leads 
us to the conclusion that the oscillations around one that 
we observe for i > 8 are numerical relics and that the true 
value of the Pi is either one or something smaller than 
one. This conclusion is based on the fact that in MPS 
simulations the transfer matrix eigenvalues that converge 
first are the dominant eigenvalues (i.e. the ones with the 
largest absolute value) so we can assume that the values 
obtained for Pi<% are much more accurate than the other 
ones. 

Thus the worst case for our purpose is when all Pi that 
are not clearly smaller than one, are actually equal to one. 
Let us investigate what we would get for the norm in this 
case. If Pi = 1, the contribution of these eigenvalues to 



Appendix C: Comparison to other PBC MPS 
algorithms 

In this appendix we will show that the algorithm |18j 
that we used to obtain all results in this work is perform- 
ing better than other recently presented approaches. 

For the sake of completeness let us first recapitulate 
the result of the comparison to the algorithm presented 
in (2D]. We have already shown in [TB] that our PBC algo- 
rithm yields a better precision. Apart from several other 
differences in these two approaches, the crucial point is 
that we allow for a variable dimension n of the dom- 
inant subspace used to approximate big powers of the 
transfer matrix. Even though this contributes a factor 
n 2 to the overall computational cost 0(n 2 D 3 ), we have 
shown in [TS] that there is no way to get rid of the fac- 
tor n if one wants to reproduce the correlation function 
throughout the entire PBC chain faithfully. If the same 
n-scanning strategy would be employed in [5D] , probably 
the same precision level could be achieved, however the 
computational cost in that algorithm would then scale 
like 0(NnD 3 ). There is an additional factor N in that 
scaling because that approach is not translationally in- 
variant. The power of n is reduced by one due to the fact 
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that the energy is minimized directly and not using the 
gradient. 

Next we would like to compare our PBC algorithm to 
the one presented in [SI] . In that work the authors sim- 
ulate the critical Quantum Ising Model by using Time 
Evolving Block Decimation to locally update a trans- 
lationally invariant MPS which is then plugged into a 
chain with PBC geometry in order to compute the energy. 
The weakness of that algorithm is that the local update of 
the MPS tensors does not take into account the boundary 
conditions: the fixed point MPS is exactly the same like 
the one obtained when trying to approximate the ground 
state of an infinite chain. In spite of this, the ground 
state energy can be approximated quite well since the 
scaling of the computational cost is only 0(nD 3 ) which 
allows the use of very large D. Unfortunately there are no 
explicit plots of the precision of the ground state energy 
in [5T] as a function of D. From the abstract and footnote 
4 of that work we deduce that the simulation that yields 
the error w 2.0 x 10~ 10 for the critical Quantum Ising 
PBC chain with N = 4800 was done with a MPS with 
bond dimension D = 200. We reach the same precision 
with D as small as 64 as can be seen in figure [l] Due to 
the higher computational cost of our algorithm D = 200 
is out of reach for us. Nonetheless we have computed an 
approximation of the ground state of the infinite chain 
with a translationally invariant MPS with D = 200 (de- 
tails of this are given below) and then plugged this MPS 
into a PBC chain geometry with N = 4800. The rela- 
tive precision that we obtained using this strategy was 
A rel E (N = 4800, D = 200) 1.39 • 10~ 10 which is in 
perfect agreement with the claim made in |21j . However, 
if we take into account the fact that a PBC simulation 
with N = 4800 and D = 200 is well in the FSES regime 
due to N < £(D = 200) w 1.910 5 , it becomes immedi- 



ately clear that with D = 200 one can in principle reach 
a much better precision than 10~ 10 . In other words, 
the results obtained in [5T] correspond to the cyan (light) 
lines in the left plot of figure [7j While this is perfectly 
fine for simulations in the FESS regime, if one is in the 
FSES regime, there is room for one or more orders of 
magnitude of improvement of the relative precision. 

There is another point worth mentioning regarding our 
PBC algorithm [TH]. In order to check the claims made 
in [21] , we needed to first approximate the ground state of 
the infinite chain with an MPS with D — 200. For this we 
used a new method called Time-Dependent Variational 
Principle j^SJ. We did this because TDVP converges 
much faster than Imaginary Time Evolution based on 
Matrix Product Operators gS] or iTEBD [37]. The rela- 
tive precision we achieved with TDVP was A re iEo(N — 
oo,D = 200) « 7.7 ■ 10~ u . We knew that this cannot 
be the best precision that can be reached with D = 200 
since in [29] we get roughly the same precision with D as 
small as 128. Thus we ran the PBC algorithm for a huge 
chain with N = 10 6 sites on top of the MPS obtained 
by TDVP. Choosing as the parameters of that algorithm 
m = 1000 and n = 100 we managed to reduce the relative 
precision to m 1.3 ■ 10 which is in perfect agreement 
with the polynomial scaling shown in figure 1 of [29j . 
The lesson learned from this approach is that TDVP re- 
sults can be improved using our PBC algorithm well in 
the FESS regime. We emphasize that running the PBC 
algorithm with small n did not yield any improvement 
to the TDVP result. Only with n as large as 100 we 
obtained the improved precision. This is quite strange 
as when we compute the energy density for the infinite 
chain, only one dominant eigenvector is used, i.e. n = 1. 
So it seems that even if additional dominant eigenvectors 
do not enter the final computation of the ground state en- 
ergy, they do have an effect during the local optimization 
procedure of the translationally invariant MPS. 



